A Multi-Rotor Drone Micro-Motion Parameter Estimation Method Based on CVMD and SVD

: It is of great signiﬁcance to detect drones in airspace due to the substantial increase and regrettable misuse in the consumer market. In this paper, we establish a micro-motion theoretical model of a drone and analyze the micro-Doppler signature of rotor targets and the ﬂicker mechanisms of the multi-rotor targets. Hence, for the target recognition problem of multi-rotor drones, a multi-rotor target micro-Doppler parameter estimation method is proposed. Firstly, a signal frequency domain segmentation method is proposed based on the complex variational mode decomposition (CVMD) to separate the high-frequency part of the high-frequency ﬂicker in the frequency domain. Secondly, for the signal after frequency domain segmentation, a ﬂicker time domain position method based on singular value decomposition (SVD) is proposed. Finally, by integrating CVMD frequency domain segmentation and SVD time domain positioning, the reconstruction of multi-rotor target scintillation at different speeds is realized, and the micro-motion parameters of rotor blades are successfully estimated. The simulation results show that the method has high accuracy in estimating the micro-motion parameters of a multi-rotor, which makes up for the shortage of the traditional method in estimating the micro-motion parameters of the multi-rotor target.


Introduction
With the development of drone technology, commercial drones have gained huge popularity in recent times, due to their ease of operation, robust performances, and low costs. Drones have been widely used in agricultural operations, video shooting, logistic transport, environmental monitoring, and rescue assistance [1]. Along with their advantages, the usability of drones also brings significant potential safety hazards, especially "black flight" and illegal activity. Small drones usually have the characteristics of low-altitude flights, small radar cross-sections (RCSs), and slow speeds ("low, small, slow"), and have very similar radar characteristics to some large birds and other jamming targets. It is difficult to distinguish these targets using traditional radar detection methods [2].
Research on the micro-Doppler effect provides a practical technical approach for identifying drones. The micro-Doppler theory was first proposed by V. C. Chen from the US Naval Laboratory [3]; it is used to represent the Doppler effect caused by the micromotion of the various components within a target. Nowadays, micro-Doppler signatures are widely used in radar target recognition and micro-motion parameter estimation [4][5][6]. The micro-motion of the drone is generated by the propeller blade rotation, and the bird is due to the flapping of the wings. These two micro-motion forms are entirely different, so there is apparent diversity in their micro-Doppler characteristics [7]. The study of micro-Doppler signatures is hence of great importance in this context [8].
The rotation of the rotor target will produce periodic modulation of the radar signal. When the rotor is at its maximum RCS position, the echo signal peaks, causing periodic flickering of the echo signal received by the radar. The flickering signal is The structure of this paper is as follows. In Section 2, the flickering mechanism of the multi-rotor target is analyzed through the micro-motion model of the rotor target. Section 3 presents a method for estimating micro-motion parameters of multi-rotor targets. Section 4 analyzes and tests the proposed algorithm. Section 5 analyzes and discusses the algorithm's performance. The last, Section 6 concludes the paper and provides an outlook for future research.

Analysis of Flickering Mechanism of Multi-Rotor Drone
In this section, the relationship between the flicker mechanism of the multi-rotor drone and the state of the rotor blades is analyzed in detail. The relationship between the motion dynamics and radar echo signal induced by the rotating blades of the drone is presented firstly. After that, the flicker mechanism of the multi-rotor micro-Doppler signature is analyzed in detail using the theoretical simulation results of the echo model. Moreover, the electromagnetic simulation is used to verify the correctness of the theoretical simulation results.

Echo Mathematical Model of a Flying Drone
The dynamic geometric relationship between the rotor drone with its rotating blades and the radar is shown in Figure 1. The radar is located at the origin of the XYZ coordinate system, R 0 is the distance from the radar center along the radar line of sight (LOS) to the drone center p 0 at the initial moment, and β is the angle between R 0 and the XOY plane. The distance between the rotor center p and the drone center p 0 is R 0p , and the angle between R 0 and R p is α. It is assumed that the drone's main body and micro-motion of its blades can be modeled using the constant radial velocity v and constant angular velocity ω p . Referring to [23], the echo signal received from the entire drone can be expressed in the time domain as s(t) =Ā · s T (t) +B · s T (t) · s M (t), in which f 0 is the radar carrier frequency and λ is the wavelength. K is the number of body scatters, A k is the scattering intensity of the kth scatterer, and R p is the distances between scatters k, p, and p 0 . α k is the angle between Op 0 and p 0 k. P is the number of drone rotor blades, Q is the number of scattering points on the blade, B pq is the scattering intensity of the scattering point q on the rotor p, and R pq is the distance between the blade scattering point q on the rotor p and the rotor axis of p. ϕ 0pq is the initial phase of the scattering point q on the rotor p, ϕ 0p is the initial phase of the rotor p, N is the number of rotor blades, and ϕ 0pq satisfies ϕ 0pq = ϕ 0p + 2π(n − 1)/N, n = 1, · · · , N − 1.
The frequency domain analysis can be obtained by taking the Fourier transform (FT) of the time domain echo: where * is the convolution operation, and δ(·) is the Dirac function. Expanding s M (t) with a Fourier series and referring to the Bessel function of the first kind, Equation (7) can be approximately expressed as where where J m (·) is the mth-order Bessel function of the first kind. Referring to Equation (9), it can be known that the Doppler frequency range of the echo signal is Figure 1. The geometric relationship between the radar and drone.

Theoretical Simulation Analysis of Flicker Mechanism and Electromagnetic Simulation Verification
In this paper, the simulation of the quad-rotor echo is based on the continuous wave (CW) radar with a carrier frequency of 10 GHz and a pulse repetition frequency of 20 kHz, and the blade scatters are regarded as a scattering line segment composed of a series of scattering points. Assume that the radar and the drone are on the same horizontal altitude; that is, the LOS angle β is 0 • , the initial distance R 0 is 50 m, and the length of the blade is 0.1 √ 2 m. In addition, to analyze the flicker features, it is assumed that the Doppler frequency shift caused by the translation motion of the drone is compensated. The simulation includes two typical motion states, hovering and horizontal uniform motion. In other drone states, considering the short observation time, the drone's blades can be approximately regarded as rotating at a uniform speed. At this time, the state of the drone can be approximately regarded as one of these two basic states. The setting parameters of the drone are shown in Table 1.
The parameter settings in Table 1 are based on the consideration of the essential force balance of the drone. (1) The sum of rotational inertia generated by the rotation of all the rotors of the drone is zero. The general control model sets the steering of the angular rotors to be the same and the steering of the adjacent rotors to be opposite. (2) The lift generated by the drone under various attitudes can balance its gravity. Assume that the quadratic sum of the rotation speed required for the hovering state of the drone is Ω. Then the blade rotating velocity satisfies the following relationship: The time-frequency spectrum obtained by the short-time Fourier transform (STFT) of the quad-rotor blades echo simulation results in Table 1 is shown in Figure 2. The frequency spectrum simulation of a single-rotor blade is mainly used for the comparative experiment of multi-rotor simulation, which consists of a group of periodic flickers (shown in Figure 2a). Comparing Figure 2a-d, it can be seen that the quad-rotor echo frequency spectrum can be seen as composed of two groups of flickers. The diagonal rotor-corresponding flickers are completely coincident and caused by the same rotating velocity and initial phase. To analyze the relationship between the relative position changes of the two groups and the initial phase of the blades, the two groups of the rotor blade's initial phase and the positions ϕ 01 , ϕ 02 , perpendicular to the radar line of sight, are introduced, as well as the evaluation parameter ξ.
where · means to round down, and this simulation ε = 90.  When ξ = 0, the two groups of flickers coincide and are consistent with the flicker of a single-rotor. With the increase of ξ, the two groups of flickers will gradually separate. Until the interval between the two groups of flickers reaches the maximum ξ = 90, all flickers are distributed at equal time intervals. The flicker interval is half of the single-rotor flicker interval. When the initial positions of the diagonal rotors are inconsistent, there will be four groups of flickers in the time spectrum (as shown in Figure 2e), which is very likely to cause severe aliasing due to the excessive flicker. The severity of aliasing depends on the difference in the initial position of the blade and the resolution of the time-frequency analysis method. When the aliasing is too severe, there will be no value in the parameter estimation.
When the drone moves at a horizontal uniform motion (equivalent to the constant radial velocity translational motion in this geometric model), there must be a difference in the rotational angular speed between the rotors. Comparing Figure 2f-h, when the initial position of the diagonal rotor blades satisfies ξ = 0, the time spectrum will be composed of two groups of flickers with different Doppler bandwidths and periods. When the initial phases of the quad-rotors are all random, the flicker will consist of two groups of flickers with varying bandwidths of Doppler in Figure 2d, in which case, the probability of severe aliasing of the flickers will be significantly increased. Moreover, we found two significant difficulties in estimating the micro-motion parameters of multi-rotor targets by comprehensively analyzing, including flicker multi-component aliasing and different component Doppler bandwidth estimations.
To study the accuracy of the theoretical simulation analysis results and expand the data sources, the electromagnetic simulation software FEKO was used to simulate the flicker features of the drone. References [15,27] built a simplified 3D model of the drone as shown in Figure 3. The simulation was divided into a fuselage component multi-attitude simulation and blade multi-attitude dynamic simulation. The blade dynamic simulation was completed by compiling and modifying the parameters of the EDITFEKO compiled file to control the rotation of the model. Finally, the time domain echo of the complete drone was constructed by the sum of the far-field electric field values of the fuselage and the blades obtained by the simulation. The simulation also adopts the state of the drone set in Table 1. Figure 4 shows the change of the time domain echo time-frequency spectrum of the single rotor with the attitude of the drone and Figure 5 offers the time-frequency spectrum of the quad-rotor in different states. Compared to Figure 2, it can be seen that the essential features of the electromagnetic simulation and the theoretical simulation are entirely consistent, which proves the reliability of the previous theoretical simulation analysis results.

Micro-Motion Parameter Estimation Algorithm Based on CVMD and SVD
For the micro-Doppler signature of multi-rotor drones, the flicker corresponding to each rotor is very likely to have aliasing. It is unrealistic to extract micro-motion parameters from the flicker with severe aliasing. The proposed algorithm takes the two groups of flickers with different frequencies, shown in Figure 5c and d; the basic process of parameter extraction is demonstrated in Figure 6. When the flickering of the quad-rotors can be distinguished, it is only necessary to make a distinction after the high and low-frequency flickers are completed (not considered in this study). In addition, the target translational motion will cause the echo signal to produce a Doppler frequency shift. At present, the mainstream method uses translation compensation to remove the influence of the target translation [19,28,29]. The algorithm in this paper is aimed at the signal that has undergone translation compensation.

CVMD Frequency Domain Segmentation
VMD was first proposed in [30], a solution process for the variational problem based on the Wiener filter, Hilbert transform, and mixed frequencies. It is assumed that the target signal consists of a specified number of band-limited intrinsic mode functions (BLIMFs) u k , each of which corresponds to a center frequency ω k and is closely clustered around the respective center frequency. Therefore, the key to realizing VMD is solving each mode's center frequency and bandwidth. The calculation result is affected by the number of decomposition layers K. The basic idea of this method is to transform the solution process of the modal bandwidth and central frequency into a constrained optimization process. The original signal is finally transformed into K BLIMFs with special sparsity and independence by iteratively searching for the optimal solution of the constrained variational model.
In view of the characteristic that VMD decomposes the signal according to the center frequency, this paper proposes using VMD to separate the high-frequency components of the high-frequency flicker in the echo signal separately. Referring to the basic principle of VMD, the VMD solution process needs to first perform the Hilbert transform on the original signal. However, the Hilbert transform can only be applied to process and preserve the complete information of the real signal. For a complex signal whose spectrum is single-sided and asymmetric, such as radar echo, half of the spectrum structure of the reconstructed signal would be lost after the Hilbert transform. Therefore, the VMD cannot directly process the radar echo signal.
In response to the defect that VMD cannot directly handle complex signals, the authors of [31,32] extended the empirical mode decomposition (EMD) method to the complex domain and successfully realized the VMD reconstruction of complex signals. The basic process of decomposing and reconstructing radar echo signals concerning CVMD is shown in Figure 7. The number of decomposed modes K needs to be modified according to different target signals.

Radar Signal
Real Signal Imaginary Signal

VMD Reconstruction
VMD Reconstruction Using the above CVMD reconstruction process, the signal is decomposed into K BLIMFs {u 1 , u 2 , · · · , u K }, and each BLIMF will correspond to a different center frequency. At this time, the spectrum obtained by FT or the time-frequency spectrum obtained by STFT can be used to sort the BLIMF according to the size of the center frequency. In this paper, the time-frequency spectrum is used to analyze the center frequency to facilitate the selection of suitable BLIMFs for combination. Finally, the frequency domain segmentation of the signal can be completed by combining the BLIMFs according to the desired purpose.

SVD Time Domain Positioning
SVD is a well-known technical means. It can be used to reduce the cross-terms in the high-resolution time-frequency distributions. It can also use the obtained singular vector as the target feature to reduce the dimension of the feature space. There are many applications in classifying "low, small, slow" objects [25,33,34]. The SVD of matrix M can be expressed as Among them, Σ is a diagonal matrix, the diagonal elements correspond to the singular values of X, and U and V represent the left singular vector and the right singular vector, respectively.
Reference [34] shows that if SVD is used for a complete time-frequency spectrum, the time and frequency information of the original spectrogram will be decoupled and projected onto the left singular vector and the right singular vector. At this time, the left singular vector U will represent the frequency information of the spectrogram, and the right singular vector V will be able to represent the time information of the spectrogram.
Considering that the singular vectors obtained by SVD are arranged according to the corresponding eigenvalues from large to small, most of the information in the time spectrum is contained in the first few singular vectors. The first three singular vectors obtained using the spectral matrix SVD of Figure 2f are shown in Figure 8. The results show that the left singular vectors can roughly reflect the micro-Doppler frequency range of the flicker band, and the frequency domain ranges of the first three left singular vectors are roughly the same. The peak of the right singular vector protrusion is very close to the time domain position of each flicker in the time-frequency spectrum. The first right singular vector is the most obvious, and the second and third right singular vectors will form multiple bifurcations near the peak of the first singular vector, which is very unfavorable for the extraction of the peak time domain position. Under careful consideration, the first singular vector obtained by SVD is applied to analysis and research. In order to obtain the temporal position of the flicker directly from the right singular vector, a thresholded peak detection method is proposed. First, a Gaussian smoothing process is performed on the first right singular vector to prevent the influence of individual Burr points to the extraction results. The length of the Gaussian smoothing window should not be too large to ensure that the influence on the elements in the original vector is minimized. Then adopt the method of peak detection to detect all the peak points in the first right singular vector. Assuming that the numerical sequence of singular vectors is {x 1 , x 2 , · · · , x n }, calculate the difference between adjacent elements ∆x i When Equation (15) is satisfied, the element can be determined as the peak point. Finally, the non-target points are eliminated through the threshold, and the threshold Th is set as Among them, δ is the threshold factor, which needs to be adjusted according to the strength of the error component, generally taking 0.6. V 1 is the first right singular vector, and max(V 1 ) represents the maximum value of V 1 . The thresholding analysis method is widely used in digital image processing, and it will also be used in the subsequent estimation of flicker micro-Doppler parameters.

Process of the Algorithm
Combined with CVMD frequency domain segmentation and SVD time domain positioning, extracting the micro-Doppler parameters of multi-rotor targets will be straightforward. The specific steps of the proposed algorithm are as follows: High-frequency flicker time domain coordinate correction. Calculate the distance between the adjacent coordinate elements ∆B.
where b i , i = 1, · · · , n are all elements in set B, and n is the number of elements in the coordinate set.
where TF f k represents the frequency domain vector whose time domain coordinate is k in the time-frequency spectrum. Then do the same threshold segmentation as formula Equation (16) with the vector FS B . At this time, the threshold generally takes δ = 0.5. After segmentation, the frequency domain width of a high-frequency flicker is the estimated value of its Doppler bandwidth B BW .
B BW = max(arg(FS B > δ · max(FS B ))) − min(arg(FS B > δ · max(FS B ))). (21) 5. The nearest neighbors difference set. Delete the point a arg min(|A−b i |) in set A that is closest to the element position in set B to obtain a low-frequency flickering time domain coordinate set C.
(a) Use the same method as step 3b to process set B to complete the missing points of set C.
Use the method of step 4 to estimate the micro-Doppler parameters C BW and ∆C of low-frequency flicker. In order to avoid the influence of the highfrequency flicker on the estimation result of the low-frequency flicker Doppler bandwidth, the flicker point set C 0 , which is relatively far away from the high-frequency flicker in set C is selected to estimate the Doppler bandwidth. Calculate the point a arg min(|C−b i |) and the corresponding shortest distance min(|C − b i |), delete all points in set C that satisfy Equation (22) to obtain set C 0 .
in which len means the window length used in STFT. (c) In order to prevent periodic missing in low-frequency flicker coordinate positioning results-use to replace ζ i as the interval evaluation factor to perform the step 3b operation to make a second correction to set C. 7.
Micro-motion parameter estimation. Referring to Equations (8) and (10), the rotor blade rotational frequency f r and blade length L can be estimated where f s represent the pulse repetition frequency, ∆X and X BW represent the estimated flicker interval and flicker Doppler bandwidth of the two groups of flickers.

Analysis of Algorithm Simulation
The proposed algorithm's purpose and the theoretical parameter estimation process have been shown in the previous subsection. In this section, we will verify the algorithm's effectiveness for estimating micro-motion parameters of multi-rotor targets. In addition, to demonstrate the specific signal processing process by the proposed algorithm, the processing process of each step of the algorithm is analyzed and displayed by taking the electromagnetic simulation results; see Figure 5c as an example.
The proposed algorithm firstly performs CVMD frequency domain segmentation on the signal after translation compensation. The key to this process is to select the appropriate number of decomposition layers K to decompose and reconstruct the original signal. In order to improve efficiency, the minimum number of decomposition layers that can separate high-frequency components is generally selected for processing. After frequency domain segmentation, the modal functions can be combined according to the time-frequency spectrum obtained by the modal function STFT. The signal will be divided into three parts: body component, high-frequency component, and intermediate frequency component. The time-frequency spectrum of the three components is shown in Figure 10. In general, the electromagnetic scattering of the body is much stronger than that of the blades, and the flicker feature in the right singular vector obtained by SVD is easily overwhelmed by the body component. Moreover, the more significant the difference between the scattering intensity of the body component and the rotor component, the more serious the contamination will be. Figure 11 shows the comparison of the first singular vector obtained before and after removing the body component. The results show that if the body component is not removed, the convex peak characterizing the flicker feature in the first right singular vector will be completely contaminated.
As shown in Figure 10c, the high-frequency components contain very few signal components. When the signal-to-noise ratio (SNR) is low, the signal pollution's influence on the time domain localization of the high-frequency flicker will be more obvious. If the coordinate set obtained from the high-frequency components is not preprocessed, the anti-noise ability of the entire parameter extraction algorithm will be significantly reduced. Figure 12a shows the flicker time domain coordinates of the high-frequency components acquired at an SNR of -2dB. Due to noise pollution, the flicker time domain location results have false detection points and missing points. The algorithm proposed in this paper adds step 3 to correct the obtained high-frequency flickering coordinates. After the first correction, the wrong coordinates will be removed (see Figure 12a). The second correction step is then performed to supplement the missed coordinates, and the final flicker coordinates detection result is shown in Figure 12c. After testing, the proposed correction strategy can obtain the correct flicker time domain coordinates as long as the obtained coordinates reflect the periodic flicker without continuous missing. The time domain coordinates of all low-frequency flickers except the aliasing flicker can be obtained through the nearest neighbor difference processing of the time domain coordinate sets of complete and high-frequency flickers. Figure 13 shows the effects of marking the high-frequency flicker coordinates and the low-frequency flicker coordinates in the time-frequency spectrum of the original signal. Due to the influence of aliasing flicker, the low-frequency flicker obtained after the difference set will be partially missing. At this time, it is also necessary to correct the low-frequency flicker coordinate set. Since the differentiated intermediate frequency components have a much more significant proportion than the high-frequency components in the original signal, its noise adaptability is more robust than the high-frequency flicker localization. On the premise that the coordinates of the high-frequency flicker can be accurately extracted, and the location results of the full flicker coordinates are accurate. The correction to the coordinates of the low-frequency flicker is aimed at complementing the missing coordinates of the aliased flicker.  Figure 14 shows the flicker reconstruction results obtained using the estimated micro-Doppler parameters. Observing the two aliased flickering areas originally in Figure 13, we found the double flickering areas in the original flickering area, which were both correctly reconstructed.

Results of Parameter Estimation
In addition to the most apparent flicker, the micro-Doppler characteristics of multirotor targets also have sinusoidal envelopes that represent the Doppler characteristics of the blade tip scattering points. The number of these sinusoids corresponds to the number of blades. Combined with the micro-Doppler curve of the flicker parameters, such as the range, the blade rotation period, and the time domain position of the initial flicker, the reconstruction of these sinusoidal envelopes can be realized. Using highfrequency flickering envelopes as an example, the mathematical expression f envelope can be expressed as f en = ±B BW sin(2π f r (t + b 1 ) + π/2) (26) Figure 15 shows the effect of drawing the reconstructed flicker and sinusoidal envelope of Figure 4c and d in the original time-frequency spectrum. It is observed that the reconstructed micro-Doppler curve is in good agreement with the micro-Doppler signatures in the original time-frequency spectrum. In order to compare the accuracy of the final extraction results, Table 2

The Influence of Noise on Parameter Estimation
To evaluate the robustness of the proposed algorithm, we test the algorithm's probability of correctly detecting the localization flicker at different SNRs and the accuracy of parameter estimation. The test uses the theoretical simulation results of type 6 and type 7 in Table 1 Figure 16 shows the probability statistics of complete flicker localization, high-frequency flicker localization, and the final reconstruction of all flickers in the parameter estimation process. The results show that the proposed algorithm can complete the flicker reconstruction when the signal SNR is larger than −2 dB.  Since the estimated flicker period accuracy can be very high when the flicker is correctly located, the primary source of error in the blade length estimation is the Doppler bandwidth estimation. Therefore, the MRE of the two groups of rotor blade lengths in Figure 17 are consistent with the Doppler bandwidth estimation results. With a comprehensive analysis of the CDP and MRE statistical results, the algorithm can effectively accomplish the micromotion parameter estimation when the SNR is larger than −2 dB.

Strengths and Weaknesses
As summarized in the following paragraphs, this paper's main strengths can be summed up as follows: (1) This paper analyzes the relationship between the multi-rotor micro-Doppler characteristics and the initial position of the rotor in detail. This work will provide a specific theoretical basis for researching the multi-rotor target micro-Doppler characteristics.
(2) The proposed CVMD frequency domain segmentation method broadens the application of VMD in radar signal processing. According to the basic idea of the proposed method, the signal with obvious frequency domain distinction can be divided into frequency domains according to different purposes. In addition, the method proposed in this paper has some shortcomings, mainly the following: (1) The frequency domain segmentation idea of CVMD relies on the stable frequency domain distinction of the signal. When the position of the target signal in the frequency domain changes according to time, the frequency domain segmentation effect of this method will decrease rapidly. (2) Generally, the rotor speed difference of the multi-rotor drone is not particularly large, and the extracted high-frequency components will contain very few signal components. Therefore, the signal-to-noise ratio that the proposed method can adapt to is limited. (3) Although the proposed method solves the flicker reconstruction in the case of weak aliasing, the aliasing situation will be dire when there are many flicker components. Therefore, it was necessary for us to introduce a time-frequency analysis method with higher resolution (especially temporal resolution) in the following research.

Conclusions
For the multi-rotor target, this paper uses the scattering point integral model to complete the construction of the theoretical echo mathematical model of the rotor drone. It completes the analysis and verification of the multi-rotor drone flicker characteristics through theoretical and electromagnetic simulations. Regarding the difficulties brought by the multi-rotor and multi-speed multi-rotor drone in estimating the micro-motion parameters, this paper takes a non-parametric parameter estimation method based on CVMD and SVD. The simulation results show that the proposed algorithm can successfully extract the micro-motion parameters of the multi-rotor target. Finally, to test the performance of the proposed algorithm, the probability of correctly detecting the position of the flicker in the time domain and the error of the parameter estimation under different SNRs were statistically analyzed. Tests show that the algorithm can complete the stable estimation of micro-motion parameters when the SNR is larger than −2dB. In follow-up research, we will consider studying the use of Doppler features of drone targets to classify their motion states. At this point, we can obtain the movement trend of the target to evaluate its risk factor. Thus, effective response strategies can be formulated in advance to avoid accidents with potential threat targets. This work will provide technical support for preventing and controlling rotary-wing drones.