In recent years, the rapid development of low-altitude economy has highlighted the importance of effective airspace management, particularly for the detection, classification, and identification of multirotor drones. Due to their small size, slow speed, and use of new materials, multirotor drones pose significant challenges for detection in urban environments. The micro-Doppler effect, resulting from micromotions of detected targets, provides a powerful means to extract motion and structural features. This paper investigates the micromotion characteristics of multirotor drones under varying wind fields, specifically focusing on gust conditions. We propose a parameter estimation method based on the segmented FFT phase difference algorithm with Hanning window to extract key micromotion features, including swing attenuation factor, amplitude, frequency, and initial phase. This approach offers a novel perspective for drone detection and identification.
The micro-Doppler effect arises from periodic movements such as rotor rotations or body oscillations in multirotor drones. When exposed to gust winds, the drone’s flight control system induces damped oscillatory responses, leading to exponential decay sinusoidal micro-Doppler signatures. These signatures serve as distinctive features for discriminating multirotor drones from other targets like birds or vehicles. In this study, we model the wind field using turbulence, gust, and wind shear models, and analyze the velocity perturbations of multirotor drones under these conditions. The results indicate that gusts cause the most significant velocity disturbances, which can be modeled as exponential decay sinusoidal signals. After effective clutter suppression, the micro-Doppler frequency of the drone’s echo signal can be expressed as:
$$ f_d(t) = \frac{2f}{c} \left[ A_0 e^{-\alpha t} \cos(2\pi f_0 t + \theta_0) + \sum_i \sin(2\pi f_i t + \theta_i) \right] $$
where $f$ is the carrier frequency, $c$ is the speed of light, $A_0$ is the amplitude, $\alpha$ is the attenuation factor, $f_0$ is the oscillation frequency, and $\theta_0$ is the initial phase. The exponential decay term represents the micromotion feature of the multirotor drone under gust influence.
To accurately estimate these parameters, we develop a segmented FFT phase difference algorithm with Hanning window. The Hanning window is chosen for its narrow main lobe and high frequency resolution, which effectively suppresses spectral leakage. The observed signal after clutter cancellation is modeled as:
$$ s(t) = A_0 e^{-\alpha t} \cos(2\pi f_0 t + \theta_0) $$
In additive white Gaussian noise environment, the sampled signal is:
$$ x(n) = s(n) + u(n), \quad n = 0, 1, \ldots, N-1 $$
where $u(n)$ is Gaussian white noise with zero mean and variance $\sigma^2$. The exponential decay sinusoidal sequence $s(n)$ is:
$$ s(n) = A_0 e^{-\alpha n / f_s} \cos \left( \frac{2\pi f_0 n}{f_s} + \theta_0 \right) $$
The signal-to-noise ratio (SNR) is defined as:
$$ \text{SNR} = \frac{\frac{1}{N} \sum_{n=0}^{N-1} s(n)^2}{\sigma^2} $$
The signal is divided into two segments. The first segment $s_1(n)$ and the second segment $s_2(n)$ are given by:
$$ s_1(n) = \frac{1}{2} A_0 e^{-\alpha n / f_s} \left[ e^{j(2\pi f_0 n / f_s + \theta_0)} + e^{-j(2\pi f_0 n / f_s + \theta_0)} \right], \quad n = 0, 1, \ldots, \frac{N}{2}-1 $$
$$ s_2(n) = s_1 \left( n + \frac{N}{2} \right) = \frac{1}{2} A_0 e^{-\alpha n / f_s} e^{-\alpha N / (2f_s)} \left[ e^{j(2\pi f_0 n / f_s + \theta_0)} e^{j\pi f_0 T} + e^{-j(2\pi f_0 n / f_s + \theta_0)} e^{-j\pi f_0 T} \right] $$
Applying the Hanning window $w(n) = \frac{1}{2} – \frac{1}{2} \cos \left( \frac{4\pi n}{N} \right)$, the windowed first segment becomes:
$$ s_{1w}(n) = \frac{1}{2} A_0 e^{-\alpha n / f_s} e^{j(2\pi f_0 n / f_s + \theta_0)} \left( 1 – \cos \left( \frac{4\pi n}{N} \right) \right) $$
After performing $\frac{N}{2}$-point FFT on $s_{1w}(n)$, the result is:
$$ S_{1w}(k) = \frac{1}{4} A_0 e^{j\theta_0} e^{\alpha(2 – T f_s) / (2f_s)} \left[ P(k) – \frac{1}{2} \left( P(k+1) + P(k-1) \right) \right] $$
where
$$ P(k) = \frac{e^{\alpha T / 2} – e^{j\pi (f_0 T – 2k)}}{e^{\alpha / f_s} – e^{j \frac{2\pi}{N} (f_0 T – 2k)}} $$
The magnitude and phase of $P(k)$ are:
$$ |P(k)| = \sqrt{ \frac{ \cosh(\alpha T / 2) – \cos(\pi (f_0 T – 2k)) }{ \cosh(\alpha / N) – \cos(\frac{2\pi}{N} (f_0 T – 2k)) } } $$
$$ \angle P(k) = \tan^{-1} \left( \frac{ \sin(\pi (f_0 T – 2k)) }{ e^{\alpha T / 2} – \cos(\pi (f_0 T – 2k)) } \right) – \tan^{-1} \left( \frac{ \sin(\frac{2\pi}{N} (f_0 T – 2k)) }{ e^{\alpha / N} – \cos(\frac{2\pi}{N} (f_0 T – 2k)) } \right) $$
Thus, $S_{1w}(k)$ can be expressed in exponential form as $S_{1w}(k) = A_{1k} e^{j\varphi_{1k}}$, with:
$$ A_{1k} = \frac{1}{4} A_0 e^{\alpha(2 – T f_s) / (2f_s)} \sqrt{ X_{re}(k)^2 + X_{im}(k)^2 } $$
$$ \varphi_{1k} = -\frac{\theta_0}{2} + \tan^{-1} \frac{X_{im}(k)}{X_{re}(k)} $$
where
$$ X_{re}(k) = |P(k)| \cos \angle P(k) – \frac{1}{2} |J(k)| \cos \angle J(k) $$
$$ X_{im}(k) = |P(k)| \sin \angle P(k) – \frac{1}{2} |J(k)| \sin \angle J(k) $$
and $J(k) = |J(k)| e^{j \angle J(k)} = \sqrt{ J_{re}^2 + J_{im}^2 } e^{j \tan^{-1} (J_{im} / J_{re}) }$, with:
$$ J_{re} = |P(k+1)| \cos \angle P(k+1) + |P(k-1)| \cos \angle P(k-1) $$
$$ J_{im} = |P(k+1)| \sin \angle P(k+1) + |P(k-1)| \sin \angle P(k-1) $$
Let $m = \lfloor f_0 T / 2 \rfloor$ be the index of the maximum magnitude line in the FFT spectrum. The phase difference between the two segments is:
$$ \Delta \varphi = \varphi_{2m} – \varphi_{1m} = \pi f_0 T – 2m\pi $$
The frequency deviation is estimated as:
$$ f_\delta = \frac{\Delta \varphi}{2\pi} \Delta f = \frac{\Delta \varphi}{\pi T} $$
where $\Delta f = 2f_s / N$ is the frequency resolution. The true frequency is then:
$$ f_0 = f_m + f_\delta = (m + \delta) \Delta f $$
with $\delta = f_0 T / 2 – m$. The attenuation factor is estimated from the ratio of the maximum magnitudes:
$$ \hat{\alpha} = \frac{2}{T} \ln \left( \frac{\hat{A}_{1m}}{\hat{A}_{2m}} \right) $$
The amplitude and initial phase are estimated using the first segment’s FFT results:
$$ \hat{A}_0 = 4 \hat{A}_{1m} e^{\frac{2f_s}{\hat{\alpha}(2 – T f_s)}} \frac{1}{ \sqrt{ \left[ |\hat{P}(m)| \cos \angle \hat{P}(m) – \frac{1}{2} |J(m)| \cos \angle J(m) \right]^2 + \left[ |\hat{P}(m)| \sin \angle \hat{P}(m) – \frac{1}{2} |J(m)| \sin \angle J(m) \right]^2 } } $$
$$ \hat{\theta}_0 = -2 \left( \hat{\varphi}_{1m} – \tan^{-1} \frac{ |\hat{P}(m)| \sin \angle \hat{P}(m) – |J(m)| \sin \angle J(m) }{ |\hat{P}(m)| \cos \angle \hat{P}(m) – |J(m)| \cos \angle J(m) } \right) $$
To validate the proposed method, we conduct simulations with parameters: $A_0 = 8$, $\alpha = 1$, $\theta_0 = \pi/3$, $f_s = 512$ Hz, $f_0 = 2$ Hz, $N = 1024$, and SNR levels of 0, 5, and 10 dB. The frequency deviation $\delta$ varies from -0.5 to 0.5. The estimation performance for each parameter is evaluated through Monte Carlo simulations with 10,000 trials.
| SNR (dB) | Method | Mean Estimated $\alpha$ | RMSE |
|---|---|---|---|
| 0 | With Hanning Window | 1.02 | 0.05 |
| Without Window | 1.15 | 0.18 | |
| 5 | With Hanning Window | 1.01 | 0.03 |
| Without Window | 1.14 | 0.17 | |
| 10 | With Hanning Window | 1.00 | 0.01 |
| Without Window | 1.13 | 0.16 |
| SNR (dB) | Method | Mean Estimated $A_0$ | RMSE |
|---|---|---|---|
| 0 | With Hanning Window | 8.05 | 0.12 |
| Without Window | 8.45 | 0.52 | |
| 5 | With Hanning Window | 8.02 | 0.08 |
| Without Window | 8.42 | 0.48 | |
| 10 | With Hanning Window | 8.00 | 0.03 |
| Without Window | 8.40 | 0.45 |
| SNR (dB) | Method | Mean Estimated $f_0$ (Hz) | RMSE (Hz) |
|---|---|---|---|
| 0 | With Hanning Window | 2.01 | 0.08 |
| Without Window | 2.05 | 0.15 | |
| 5 | With Hanning Window | 2.00 | 0.05 |
| Without Window | 2.04 | 0.12 | |
| 10 | With Hanning Window | 2.00 | 0.02 |
| Without Window | 2.03 | 0.10 |
| SNR (dB) | Method | Mean Estimated $\theta_0$ (rad) | RMSE (rad) |
|---|---|---|---|
| 0 | With Hanning Window | 1.05 | 0.10 |
| Without Window | 1.12 | 0.25 | |
| 5 | With Hanning Window | 1.03 | 0.07 |
| Without Window | 1.10 | 0.22 | |
| 10 | With Hanning Window | 1.01 | 0.04 |
| Without Window | 1.08 | 0.18 |
The simulation results demonstrate that the proposed method with Hanning window significantly improves estimation accuracy for all parameters, especially the attenuation factor and amplitude, which are key micromotion features of multirotor drones under gust conditions. The Hanning window effectively suppresses spectral leakage, leading to lower RMSE values across all SNR levels. For frequency estimation, the windowed method outperforms the non-windowed approach when the frequency deviation $\delta$ is within $[-0.3, 0.4]$. At larger deviations, noise effects may cause errors in identifying the maximum magnitude line, but overall performance remains satisfactory. The initial phase estimation also benefits from the windowing technique, with improved accuracy and reduced sensitivity to frequency deviations.
In conclusion, this study establishes that multirotor drones under gust winds exhibit exponential decay sinusoidal micro-Doppler signatures, which can serve as effective micromotion features for detection and identification. The segmented FFT phase difference algorithm with Hanning window provides high-precision parameter estimation, enabling reliable discrimination of multirotor drones from other targets. Future work will focus on experimental validation using real-world data and extending the method to multi-drone scenarios.